﻿using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Linq;
using System.Text;
using System.Windows.Forms;
using DotSpatial.Controls;
using DotSpatial.Topology;
using DotSpatial.Symbology;
using DotSpatial.Data;
using System.Collections;
using DotSpatial.Projections;
using System.IO;

namespace LiDARAspect
{
    public partial class LiDARAspect : Form
    {
        public DotSpatial.Data.Raster DEMRaster = new DotSpatial.Data.Raster();
        public DotSpatial.Topology.KDTree.KdTree _tree = new DotSpatial.Topology.KDTree.KdTree(2);//define a tree
        public DotSpatial.Data.FeatureSet LiDARFeatureSet = new DotSpatial.Data.FeatureSet(FeatureType.Point);//LiDAR points or Synthetic Points
        public DotSpatial.Data.Raster GridFile = new DotSpatial.Data.Raster();
        public DotSpatial.Data.FeatureSet OutLet = new DotSpatial.Data.FeatureSet();
        public FeatureSet SyntheticPoints = new FeatureSet(FeatureType.Point);
        public FeatureSet SlopePoint = new FeatureSet(FeatureType.Point);
        public FeatureSet AspectPoint = new FeatureSet(FeatureType.Point);
        public FeatureSet LASPoint = new FeatureSet(FeatureType.Point);
        public FeatureSet NetLine = new FeatureSet(FeatureType.Line);
        public int numofshape;
        public int nNodes = 100000;
        public Node[] Nodes;
        public ArrayList randonList = new ArrayList();
        public LiDARAspect()
        {
            InitializeComponent();
            map1.GeoMouseMove += Map_GeoMouseMove;
        }

        private void Map_GeoMouseMove(object sender, GeoMouseArgs e)
        {
            tssCoordinateLabel.Text = String.Format("X: {0}, Y: {1}", e.GeographicLocation.X, e.GeographicLocation.Y);
        }

        private void loadDEMToolStripMenuItem_Click(object sender, EventArgs e)
        {
            map1.ProjectionModeDefine = ActionMode.Never;
            IMapRasterLayer RLayer = map1.AddRasterLayer();

            //Try to not use the AddRasterLayer() function by using raster open instead.
            // map1.ProjectionModeReproject = ActionMode.Never;
            if (RLayer == null)
            {
                MessageBox.Show("Please select a raster layer.");
            }
            else
            {

               // DEMRaster = new (IRaster)RLayer.DataSet;
                IRaster r = RLayer.DataSet;
                GridFile = r as Raster;
                //DEMRaster

            }
            DEMRaster.Projection = map1.Projection;
            RLayer.Symbolizer.ShadedRelief.IsUsed = true;
            RLayer.Symbolizer.ShadedRelief.ElevationFactor = -4;
            RLayer.WriteBitmap();
        }

        private void addOutletToolStripMenuItem_Click(object sender, EventArgs e)
        {
            OpenFileDialog ofd = new OpenFileDialog();
            ofd.Filter = "SHP|*.shp";
            if (ofd.ShowDialog() != DialogResult.OK) return;
            OutLet.Open(ofd.FileName);
            map1.AddLayer(ofd.FileName);
        }

        private void creatSyntheticPointsToolStripMenuItem_Click(object sender, EventArgs e)
        {
            if (GridFile == null) return;
            Random r = new Random();
            ArrayList list = new ArrayList(nNodes);

            for (int i = 0; i < nNodes; i++)
            {
                list.Add(i);
            }

            for (int i = 0; i < nNodes; i++)
            {
                int index = r.Next(0, list.Count);
                randonList.Add(list[index]);
                list.RemoveAt(index);
            }
            int Row = (int)GridFile.Maximum;
            int Col = (int)GridFile.Minimum;

            // DotSpatial.Symbology.RandomExt rnd ;//= new RandomExt();
            nNodes += OutLet.Features.Count;

            Nodes = new Node[nNodes];

            DotSpatial.Data.RcIndex RI = new DotSpatial.Data.RcIndex();

            double[] dbl = new double[2];
            for (int i = 0; i < nNodes; i++)
            {
                Nodes[i] = new Node();
                if (i < OutLet.Features.Count)
                {
                    Nodes[i].X = OutLet.Features[i].Coordinates[0].X;
                    Nodes[i].Y = OutLet.Features.ElementAt(i).Coordinates[0].Y;
                }
                else
                {
                    Nodes[i].X = RandBetween(GridFile.Xllcenter, GridFile.Xllcenter + GridFile.NumColumns * GridFile.CellWidth, i);
                    //Nodes[i].Y = RandBetween(GridFile.Yllcenter,GridFile.Yllcenter + GridFile.NumRows*GridFile.CellHeight,i);
                    Nodes[i].Y = GridFile.Yllcenter + (GridFile.Extent.Height / (double)nNodes) * i;
                }

                dbl[0] = Nodes[i].X;
                dbl[1] = Nodes[i].Y;
                Coordinate cod = new Coordinate(dbl[0], dbl[1]);
                DotSpatial.Topology.Point pt = new DotSpatial.Topology.Point(cod);
                _tree.Insert(dbl, Nodes[i]);
                pt.X = Nodes[i].X;
                pt.Y = Nodes[i].Y;
                DotSpatial.Data.IFeature f = SyntheticPoints.AddFeature(pt);//Add the feature to the FeatureSet
                RI = GridFile.Bounds.ProjToCell(cod);
                try
                {
                    Nodes[i].Z = GridFile.Value[RI.Row, RI.Column];
                }
                catch
                {
                    if (i == 0) Nodes[i].Z = 0.0;
                    else
                    {
                        Nodes[i].Z = Nodes[i - 1].Z;
                    }
                }
            }
            SyntheticPoints.Projection = map1.Projection;
            MapPointLayer l = new MapPointLayer(SyntheticPoints);
            l.Name = "SyntheticPoints";
            map1.Layers.Add(l);
        }

        /// <summary>
        /// Add the prepared LiDAR shape file and build up the KD-Tree
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void addLiDARShpToolStripMenuItem_Click(object sender, EventArgs e)
        {
            //loadDEMToolStripMenuItem_Click( sender, e);
            addOutletToolStripMenuItem_Click(sender, e);
            OpenFileDialog ofd = new OpenFileDialog();
            ofd.Filter = "SHP|*.shp";
            if (ofd.ShowDialog() != DialogResult.OK) return;
            LiDARFeatureSet.Open(ofd.FileName);
            //map1.AddLayer(ofd.FileName);
            nNodes = LiDARFeatureSet.Features.Count;
            nNodes += OutLet.Features.Count;
            Nodes = new Node[nNodes];
            double[] dbl = new double[2];
            for (int i = 0; i < nNodes; i++)
            {
                Nodes[i] = new Node();
                if (i < OutLet.Features.Count)
                {
                    Nodes[i].X = OutLet.Features[i].Coordinates[0].X;
                    Nodes[i].Y = OutLet.Features.ElementAt(i).Coordinates[0].Y;
                    //Nodes[i].Z = 0;
                }
                else
                {
                    Nodes[i].X = LiDARFeatureSet.Features[i - 1].Coordinates[0].X;
                    Nodes[i].Y = LiDARFeatureSet.Features[i - 1].Coordinates[0].Y;
                    Nodes[i].Z = LiDARFeatureSet.Features[i - 1].Coordinates[0].Z;
                }

                dbl[0] = Nodes[i].X;
                dbl[1] = Nodes[i].Y;
                //dbl[2] = Nodes[i].Z;
                //dbl[3] = i;
                Coordinate cod = new Coordinate(dbl[0], dbl[1]);
                DotSpatial.Topology.Point pt = new DotSpatial.Topology.Point(cod);
                if (_tree.Search(dbl) == null)
                {
                    _tree.Insert(dbl, Nodes[i]);
                }
                    pt.X = Nodes[i].X;
                    pt.Y = Nodes[i].Y;
                    pt.Z = Nodes[i].Z;
                    DotSpatial.Data.IFeature f = SyntheticPoints.AddFeature(pt);//Add the feature to the FeatureSet
                //}
                //catch
                //{
                //    MessageBox.Show( i.ToString(),Nodes[i].X.ToString());
                //    //Nodes[i].Y.ToString()Nodes[i].Z.ToString());
                //}
            }
            SyntheticPoints.Projection = map1.Projection;
            MapPointLayer l = new MapPointLayer(SyntheticPoints);
            l.Name = "SyntheticPoints";
            map1.Layers.Add(l);
        }
        private double RandBetween(double Low, double High, int i)
        {
            double randomF = (Convert.ToDouble(randonList[i - 1].ToString())) / (double)(randonList.Count);
            return (High - Low) * randomF + Low;
        }
        /// <summary>
        /// This function used a recursive method to determine whether a node has a downnode or not
        /// </summary>
        /// <param name="CurrentNode"></param>
        /// <param name="K"></param>
        private void ProcessSlope(Node CurrentNode, int K)
        {
            double[] d = new double[nNodes - 1];
            int[] keys = new int[nNodes - 1];
            double slope;
            double[] Location = new double[2];
            Location[0] = CurrentNode.X;
            Location[1] = CurrentNode.Y;
            //int K = 20; //the number of neighbors;
            object[] Neighbors = new object[K];
            Neighbors = _tree.Nearest(Location, K);
            for (int i = 0; i < K; i++)
            {
                Node Neighbor = new Node();
                Neighbor = (Node)(Neighbors[i]);
                if (Neighbor.Z > CurrentNode.Z)
                {
                    double oldSlope = Neighbor.Slope();
                    slope = Neighbor.GetSlope(CurrentNode);
                    if (slope > oldSlope)
                    {
                        Neighbor.DownNode = CurrentNode;
                    }
                    if (Neighbor.Visited == false) ProcessSlope(Neighbor, K);
                }
            }
            CurrentNode.Visited = true;
            DotSpatial.Topology.Point HigherPoint = new DotSpatial.Topology.Point();
            HigherPoint.X = CurrentNode.X;
            HigherPoint.Y = CurrentNode.Y;
            HigherPoint.Z = CurrentNode.Z;
            DotSpatial.Data.IFeature f = SlopePoint.AddFeature(HigherPoint);
        }

        /// <summary>
        /// Using the slope based climbing up algorithm for the stream delineation
        /// </summary>
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void slopeDelineationToolStripMenuItem_Click(object sender, EventArgs e)//Test the slope based delineation
        {

            ProcessSlope(Nodes[0], 5000);
            DotSpatial.Topology.Point p = new DotSpatial.Topology.Point();
            for (int i = 0; i < Nodes.Length; i++)
            {
                if (Nodes[i].DownNode != null)
                {
                    DotSpatial.Topology.Point[] LinePoints = new DotSpatial.Topology.Point[2];
                    p.X = Nodes[i].X;
                    p.Y = Nodes[i].Y;
                    p.Z = Nodes[i].Z;
                    LinePoints[0] = new DotSpatial.Topology.Point(Nodes[i].X, Nodes[i].Y, Nodes[i].Z);
                    LinePoints[1] = new DotSpatial.Topology.Point(Nodes[i].DownNode.X, Nodes[i].DownNode.Y, Nodes[i].DownNode.Z);
                    DotSpatial.Topology.LineString ls = new LineString(LinePoints);
                    IFeature LineFeature = NetLine.AddFeature(ls);
                }
            }
            NetLine.Projection = map1.Projection;
            MapLineLayer Line = new MapLineLayer(NetLine);
            Line.Name = "Stream Network";
            map1.Layers.Add(Line);
        }
        /// <summary>
        /// Compute the aspect of the given points
        /// </summary>
        /// should based on the kd-tree
        /// <param name="sender"></param>
        /// <param name="e"></param>
        private void aspectShowToolStripMenuItem_Click(object sender, EventArgs e)
        {
            ProcessAspect(Nodes[1], 80);//Assume that there are eight points in 8 directions.
            for (int i = 0; i < OutLet.Features.Count; i++)
            {
                ProcessAspect(Nodes[i], 8);
            }
            AspectPoint.DataTable.Columns.Add(new DataColumn("Aspect", typeof(double)));

            for (int i = 0; i < Nodes.Length-2; i++)
            {
                if (Nodes[i].DownNode != null)
                {
                    DotSpatial.Topology.Point AspectP = new DotSpatial.Topology.Point();
                    AspectP.X = Nodes[i].X;
                    AspectP.Y = Nodes[i].Y;
                    AspectP.Z = Nodes[i].Z;
                    DotSpatial.Data.IFeature f = AspectPoint.AddFeature(AspectP);
                    f.DataRow.BeginEdit();
                    f.DataRow["Aspect"] = Nodes[i].GetAspect(Nodes[i].DownNode);
                    f.DataRow.EndEdit();
                }
            }
            MapPointLayer AspectPL = new MapPointLayer(AspectPoint);
            AspectPoint.Projection = map1.Projection;
            map1.Layers.Add(AspectPL);
        }

        private void ProcessAspect(Node CurrentNode, int K)
        {
            double[] d = new double[nNodes - 1];
            int[] keys = new int[nNodes - 1];
            double Aspect;
            double[] Location = new double[2];
            Location[0] = CurrentNode.X;
            Location[1] = CurrentNode.Y;
            //int K = 20; //the number of neighbors;
            object[] Neighbors = new object[K];
            Neighbors = _tree.Nearest(Location, K);
            for (int i = 0; i < K; i++)
            {
                Node Neighbor = new Node();
                Neighbor = (Node)(Neighbors[i]);
                if (Neighbor.Z > CurrentNode.Z)
                {
                    double oldAspect = Neighbor.Aspect();
                    Aspect = Neighbor.GetSlope(CurrentNode);
                    if (Aspect > oldAspect)
                    {
                        Neighbor.DownNode = CurrentNode;
                    }
                    if (Neighbor.Visited == false) ProcessAspect(Neighbor, K);
                }
            }
            CurrentNode.Visited = true;
            DotSpatial.Topology.Point HigherPoint = new DotSpatial.Topology.Point();
            HigherPoint.X = CurrentNode.X;
            HigherPoint.Y = CurrentNode.Y;
            HigherPoint.Z = CurrentNode.Z;
            DotSpatial.Data.IFeature f = AspectPoint.AddFeature(HigherPoint);
        }

        private void loadLiDARPointToolStripMenuItem_Click(object sender, EventArgs e)
        {
            OpenFileDialog ofd = new OpenFileDialog();
            ofd.Filter = "LAS|*.las";
            if (ofd.ShowDialog() != DialogResult.OK) return;
            string filename = ofd.FileName;
            LASReader L = new LASReader(filename);
            //Loop through and read each of the points in the LAS reader
            int PTN = (int)L.header.NumberOfPointRecords;
            Nodes = new Node[PTN];
            //*****************code for create the field for store the informatiom from LiDAR Point**************//
            //LASPoint.DataTable.Columns.Add(new DotSpatial.Data.Field("Id", DotSpatial.Data.FieldDataType.String));
            //LASPoint.DataTable.Columns.Add(new DotSpatial.Data.Field("X", DotSpatial.Data.FieldDataType.Double));
            //LASPoint.DataTable.Columns.Add(new DotSpatial.Data.Field("Y", DotSpatial.Data.FieldDataType.Double));
            //LASPoint.DataTable.Columns.Add(new DotSpatial.Data.Field("Z", DotSpatial.Data.FieldDataType.Double));
            //LASPoint.DataTable.Columns.Add(new DotSpatial.Data.Field("GPSTime", DotSpatial.Data.FieldDataType.Double));
            //LASPoint.DataTable.Columns.Add(new DotSpatial.Data.Field("ScanAngle", DotSpatial.Data.FieldDataType.String));

            LASPoint.DataTable.Columns.Add(new DataColumn("Id", typeof(String)));
            LASPoint.DataTable.Columns.Add(new DataColumn("X", typeof(double)));
            LASPoint.DataTable.Columns.Add(new DataColumn("Y", typeof(double)));
            LASPoint.DataTable.Columns.Add(new DataColumn("Z", typeof(double)));
            LASPoint.DataTable.Columns.Add(new DataColumn("GPSTime", typeof(double)));
            LASPoint.DataTable.Columns.Add(new DataColumn("ScanAngle", typeof(String)));
            ////**************add the Aspect and Slope fields **********///////////
            LASPoint.DataTable.Columns.Add(new DataColumn("Slope", typeof(double)));
            LASPoint.DataTable.Columns.Add(new DataColumn("Aspect", typeof(double)));
            
            for (uint i = 0; i < PTN-1; i++)
            {
                DotSpatial.Topology.Point p = new DotSpatial.Topology.Point();
                Nodes[i] = new Node();
                Nodes[i].X = p.X = L.GetPointN(i).GetX();
                Nodes[i].Y = p.Y = L.GetPointN(i).GetY();
                Nodes[i].Z = p.Z = L.GetPointN(i).GetZ();
                DotSpatial.Data.IFeature f = LASPoint.AddFeature(p);
                //*****************code for record the information to data table**************//
                #region
                //f.DataRow.BeginEdit();
                f.DataRow["ID"] = i.ToString();
                f.DataRow["X"] = p.X;
                f.DataRow["Y"] = p.Y;
                f.DataRow["Z"] = p.Z;
                f.DataRow["GPSTime"] = L.GetPointN(i).GetGPSTime();
                f.DataRow["ScanAngle"] = L.GetPointN(i).GetScanAngle();
                f.DataRow.EndEdit();
                #endregion
                double[] dbl = new double[2];
                dbl[0] = Nodes[i].X;
                dbl[1] = Nodes[i].Y;
                try
                {
                    _tree.Insert(dbl, Nodes[i]);//build the kd-tree, LiDAR points are too dense to build up the tree.
                }
                catch
                {
                }
            }

            String Directory = Application.StartupPath;
            String SaveFilename = filename + ".shp";
            LASPoint.SaveAs(SaveFilename, true);
            LASPoint.Open(SaveFilename);
            LASPoint.Projection = map1.Projection;
            map1.AddLayer(LASPoint.Filename);

            /*
            //calculate slope and aspect for each point
            ProcessAspect(Nodes[1], 80);//Assume that there are eight points in 8 directions.
            ProcessSlope(Nodes[1], 80);
   
            AspectPoint.DataTable.Columns.Add(new DataColumn("Aspect", typeof(double)));
            AspectPoint.DataTable.Columns.Add(new DataColumn("Slope", typeof(double)));
            for (int i = 0; i < Nodes.Length - 2; i++)
            {
                if (Nodes[i].DownNode != null)
                {
                    DotSpatial.Topology.Point AspectP = new DotSpatial.Topology.Point();
                    AspectP.X = Nodes[i].X;
                    AspectP.Y = Nodes[i].Y;
                    AspectP.Z = Nodes[i].Z;
                    DotSpatial.Data.IFeature f = AspectPoint.AddFeature(AspectP);
                    f.DataRow.BeginEdit();
                    f.DataRow["Aspect"] = Nodes[i].GetAspect(Nodes[i].DownNode);
                    f.DataRow["Slope"] = Nodes[i].GetSlope(Nodes[i].DownNode);
                    f.DataRow.EndEdit();
                }
                else
                {
                    DotSpatial.Topology.Point AspectP = new DotSpatial.Topology.Point();
                    DotSpatial.Data.IFeature f = AspectPoint.AddFeature(AspectP);
                    f.DataRow.BeginEdit();
                    f.DataRow["Aspect"] = 0;
                    f.DataRow["Slope"] = 0;
                }
            }
            String AspectFilename = filename +"Aspect" + ".shp";
            AspectPoint.SaveAs(AspectFilename, true);
            AspectPoint.Open(AspectFilename);
            AspectPoint.Projection = map1.Projection;
            map1.AddLayer(AspectPoint.Filename);
            //MapPointLayer LasP = new MapPointLayer(LASPoint);
            //LASPoint.Projection = map1.Projection;
            //map1.Layers.Add(LasP);   
             */
        }

    }

    public class Node
    {
        public double X, Y, Z;
        public Node DownNode;
        public bool Visited;
        //public double Aspect;
        public double Aspect()
        {
            if (DownNode == null)
                return 0;
            else return GetAspect(DownNode);
        }
        public double Slope()
        {
            if (DownNode == null)
                return 0;
            else return GetSlope(DownNode);
        }
        public double GetSlope(Node DownNode)
        {
            double dx, dy, dist; ;
            dx = X - DownNode.X;
            dy = Y - DownNode.Y;
            dist = Math.Sqrt(dx * dx + dy * dy);
            return (Z - DownNode.Z) / dist;
        }

        /// <summary>
        /// Compute the aspect based on the angle in 360 degrees
        /// </summary>
        /// <param name="DownNode"></param>
        /// <returns></returns>
        public double GetAspect(Node DownNode)
        {
            double dx,dy,dz,dist,Aspect;
            dx = X - DownNode.X;
            dy = Y - DownNode.Y;
            dz = Z - DownNode.Z;
            dist = Math.Sqrt(dx*dx + dy*dy);
            Aspect = 0;
            if ((dx <= 0) && (dy < 0))
            {
                Aspect = Math.Asin(Math.Abs(dx) / dist) * 180 / Math.PI;
            }
            if ((dx < 0) && (dy >= 0 ))
            {
                Aspect = 90 + Math.Asin(Math.Abs(dy) / dist) * 180 / Math.PI;
            }
            if ((dx >= 0) && (dy > 0))
            {
                Aspect = 180 + Math.Asin(Math.Abs(dx) / dist) * 180 / Math.PI;
            }
            if ((dx > 0) && (dy <= 0))
            {
                Aspect = 270 + Math.Asin(Math.Abs(dy) / dist) * 180 / Math.PI;
            }
            //Aspect = (180/Math.PI) * Math.Atan((dz/dy)/(dz/dx));//this is wrong, because it is larger than 360!
            //Aspect = 180 -(360/2*PI) * (Math.Atan((dz / dy) / (dz / dx))) + 90 * (dz / dx) / Math.Abs(dz / dx);
            return Aspect;
        }
        //Other method to calculate the hydrology slope and aspect
    }
}

/*
 * need to add a new algorithm on using the dropping down method
 * One good way can be among the points that get from the climbing up point, to process a dropping down delineation
 */
